def InputSpecifier():
    import re
    ifilter1=[]
    # 0. Fractional position
    ifilter1.append(re.compile(r"pos\s+(\S+)\s+(\S+)\s+(\S+)"))
    # 1. Direction
    ifilter1.append(re.compile(r"dir\s+(\S+)\s+(\S+)\s+(\S+)"))
    # 2. Recoil energy
    ifilter1.append(re.compile(r"rec\s+(\d+)\s*(\D+)"))
    # 3. Source file name
    ifilter1.append(re.compile(r"file\s+(\S+)"))
    # 4. Atomic mass in AMU
    ifilter1.append(re.compile(r"mass\s+(\S+)"))
    for item in ifilter1:
        yield item

def iround(x):
    return int(round(x) - 0.5)+(x>0)

def readinp(fp,log_fp):
    import re
    while True:
        eatline=fp.readline()
        neat=len(eatline.strip())
        if not eatline:
            break
        log_fp.write(" "+eatline.strip()+"\n")
        for i,if1 in enumerate(InputSpecifier()):
            f1=if1.match(eatline)
            if bool(f1):
                if i == 0:
                    # Fractional position
                    px=float(f1.group(1))
                    py=float(f1.group(2))
                    pz=float(f1.group(3))
                elif i == 1:
                    # Direction
                    dx=float(f1.group(1))
                    dy=float(f1.group(2))
                    dz=float(f1.group(3))
                elif i == 2:
                    # Recoil energy
                    rec=float(f1.group(1))
                    if bool(re.search(r"[kK]",f1.group(2))):
                        rec=rec*1.0e+3
                    if bool(re.search(r"[mM]",f1.group(2))):
                        rec=rec*1.0e+6
                elif i == 3:
                    # source file
                    fn=f1.group(1).strip()
                elif i == 4:
                    mass=float(f1.group(1))
                else:
                    pass
                break
            else:
                pass
    return fn,px,py,pz,dx,dy,dz,rec,mass

def setinit():
    import sys
    import time

    if __name__=="__main__":

        try:
            ifn=sys.argv[1]
        except:
            print("Specify input file name.")
            sys.exit(2)

        try:
            ofn=sys.argv[2]
        except:
            print("Specify output file name.")
            sys.exit(2)

    else:

        ifn=input(" - Specify input file name")
        ofn=input(" - Specify output file name")

    t0=time.time()
    lfn="log%.0f" % t0

    ifp=open(ifn.strip(),'r')
    lfp=open(lfn.strip(),'w')
    ofp=open(ofn.strip(),'w')

    str1=" "+time.asctime(time.gmtime(t0))+"\n"
    str2=" Input file name: %s\n Output file name: %s\n Log file name: %s\n" % (ifn,ofn,lfn)
    sys.stdout.write(str1)
    sys.stdout.write(str2)
    lfp.write(str1)
    lfp.write(str2)
    lfp.write("\n -Input specifier..\n")
    sfn,px,py,pz,dx,dy,dz,rec,amass=readinp(ifp,lfp)
    lfp.write("\n -Find nearest atoms..")    
    sfp=open(sfn,'r')
    
    # Time step
    eatline=sfp.readline()
    eatline=sfp.readline()
    # Number of atoms
    eatline=sfp.readline()
    eatline=sfp.readline()
    natom=int(eatline)
    # Box bounds
    eatline=sfp.readline()
    eatline=sfp.readline()
    bound=eatline.split()
    bx1=float(bound[0])
    bx2=float(bound[1])
    eatline=sfp.readline()
    bound=eatline.split()
    by1=float(bound[0])
    by2=float(bound[1])
    eatline=sfp.readline()
    bound=eatline.split()
    bz1=float(bound[0])
    bz2=float(bound[1])
    # Find nearest position
    xx=bx1+(bx2-bx1)*px
    yy=by1+(by2-by1)*py
    zz=bz1+(bz2-bz1)*pz
    # Dummy line
    eatline=sfp.readline()
    dist=1.0e+6
    nearest=-1

    showprocess=0.0
    processdelta=0.1
    log_str=" -complete reading %3.2f %%\n"
    for i in range(natom):
        curprocess=float(i)/float(natom)
        if curprocess>=showprocess:
            lfp.write(log_str % (showprocess*100))
            sys.stdout.write(log_str % (showprocess*100))
            showprocess+=processdelta
        eatline=sfp.readline()
        eatword=eatline.split()
        atom_id=int(eatword[0])
        atom_type=int(eatword[1])
        x=float(eatword[2])
        y=float(eatword[3])
        z=float(eatword[4])
        cdist=(xx-x)**2+(yy-y)**2+(zz-z)**2
        if cdist < dist:
            nx=x
            ny=y
            nz=z
            nearest=atom_id
            dist=cdist
    if nearest<0:
        raise Exception("Fatal error occurs.")
    # Find vx, vy, vz
    import math
    amu2ev=931.494e+6 # eV/c2
    vc=2.99792458e+8 # m/sec
    vc=vc*0.01 # A/ps
    # rec in eV, amass in amu
    k=math.sqrt((2.0*rec*vc**2)/(amass*amu2ev*(dx**2+dy**2+dz**2)))
    vx=k*dx
    vy=k*dy
    vz=k*dz

    # Print result
    str1=" Nearest atom index=%d px=%.3f py=%.3f pz=%.3f (distance %.3f from %.3f %.3f %.3f)\n"
    str2=" E: %.3f eV, vx: %.3e, vy: %.3e, vz: %.3e in A/ps (Dir: %d %d %d)\n"
    ofp.write(str1 % (nearest,nx,ny,nz,dist,xx,yy,zz))
    ofp.write(str2 % (rec,vx,vy,vz,dx,dy,dz))
    lfp.write(str1 % (nearest,nx,ny,nz,dist,xx,yy,zz))
    lfp.write(str2 % (rec,vx,vy,vz,dx,dy,dz))
    sys.stdout.write(str1 % (nearest,nx,ny,nz,dist,xx,yy,zz))
    sys.stdout.write(str2 % (rec,vx,vy,vz,dx,dy,dz))

    sfp.close()
    lfp.close()
    ofp.close()

setinit()
